8.3 workflow全流程R包

8.3.1 Biomod2(script)

8.3.1.1 单物种建模

8.3.1.1.1 建模引用

## Citation:
## @book{
##    title={Habitat Suitability and Distribution Models: With Applications in R},
##    author={Guisan, A. and Thuiller, W. and Zimmermann, N.E.},
##    isbn={9780521758369},
##    series={Ecology, Biodiversity and Conservation},
##    year={2017},
##    publisher={Cambridge University Press}
## }

8.3.1.1.2 加载工作环境及数据

## load the required packages
library(biomod2)
library(ggplot2)
library(gridExtra)
library(raster)
library(rasterVis)
## read data ----
ProLau_occ <- read.csv('../data/ProLau_occ.csv')
summary(ProLau_occ)

bioclim_ZA_sub <- 
  raster::stack(
    c(
      bio_5  = '../data/worldclim_ZA/worldclim_ZA_bio_5.asc',
      bio_7  = '../data/worldclim_ZA/worldclim_ZA_bio_7.asc',
      bio_11 = '../data/worldclim_ZA/worldclim_ZA_bio_11.asc',
      bio_19 = '../data/worldclim_ZA/worldclim_ZA_bio_19.asc'
    )
  )
bioclim_ZA_sub
## format the data ---- ##注意这里构建缺失点数据会构建两套；#还没搞清楚为什么构建了两套；
ProLau_data <- 
  BIOMOD_FormatingData(
    resp.var = ProLau_occ['Protea.laurifolia'],
    resp.xy = ProLau_occ[, c('long', 'lat')],
    expl.var = bioclim_ZA_sub,
    resp.name = "Protea.laurifolia",
    PA.nb.rep = 2,
    PA.nb.absences = 500,
    PA.strategy = 'random'
  )
## 查看构建的数据集，为S4形式；
## formatted object summary
ProLau_data
## plot of selected pseudo-absences
plot(ProLau_data) ##包含三张图：实际分布和两张构建的缺失点分布

8.3.1.1.3 构筑多模型建模

## define individual models options ---- 
ProLau_opt <- 
  BIOMOD_ModelingOptions(
    GLM = list(type = 'quadratic', interaction.level = 1),
    GBM = list(n.trees = 1000),
    GAM = list(algo = 'GAM_mgcv')
  )

## run the individual models ----
ProLau_models <- 
  BIOMOD_Modeling(
    data = ProLau_data,
    models = c("GLM", "GBM", "RF", "GAM"),
    models.options = ProLau_opt,
    NbRunEval = 2,
    DataSplit = 80,
    VarImport = 3,
    modeling.id = "demo1"
  )

8.3.1.1.3 建模评估:AUC/TSS

## get models evaluation scores
ProLau_models_scores <- get_evaluations(ProLau_models)
## ProLau_models_scores is a 5 dimension array containing the scores of the models
dim(ProLau_models_scores)
dimnames(ProLau_models_scores)
## plot models evaluation scores
models_scores_graph(  ##多个模型的评估auc和tss
  ProLau_models, 
  by = "models", 
  metrics = c("ROC","TSS"), 
  xlim = c(0.5,1), 
  ylim = c(0.5,1)
)
models_scores_graph(  ##两次模型运行的结果评估
  ProLau_models, 
  by = "cv_run" , 
  metrics = c("ROC","TSS"), 
  xlim = c(0.5,1), 
  ylim = c(0.5,1)
)
models_scores_graph(  ##不同缺失点数据集的结果评估
  ProLau_models, 
  by = "data_set", 
  metrics = c("ROC","TSS"), 
  xlim = c(0.5,1), 
  ylim = c(0.5,1)
)
## check variable importance
(ProLau_models_var_import <- get_variables_importance(ProLau_models))
## make the mean of variable importance by algorithm ##统计汇总不同模型的TSS、AUC等；
apply(ProLau_models_var_import, c(1,2), mean)

8.3.1.1.4 建模评估:响应曲线

## individual models response plots  ##绘制不同模型运行后的响应曲线
ProLau_glm <- BIOMOD_LoadModels(ProLau_models, models='GLM')
ProLau_gbm <- BIOMOD_LoadModels(ProLau_models, models='GBM')
ProLau_rf <- BIOMOD_LoadModels(ProLau_models, models='RF')
ProLau_gam <- BIOMOD_LoadModels(ProLau_models, models='GAM')

glm_eval_strip <- 
  biomod2::response.plot2(
    models  = ProLau_glm,
    Data = get_formal_data(ProLau_models,'expl.var'), 
    show.variables= get_formal_data(ProLau_models,'expl.var.names'),
    do.bivariate = FALSE,
    fixed.var.metric = 'median',
    legend = FALSE,
    display_title = FALSE,
    data_species = get_formal_data(ProLau_models,'resp.var')
  )

gbm_eval_strip <- 
  biomod2::response.plot2(
    models  = ProLau_gbm,
    Data = get_formal_data(ProLau_models,'expl.var'), 
    show.variables= get_formal_data(ProLau_models,'expl.var.names'),
    do.bivariate = FALSE,
    fixed.var.metric = 'median',
    legend = FALSE,
    display_title = FALSE,
    data_species = get_formal_data(ProLau_models,'resp.var')
  )

rf_eval_strip <- 
  biomod2::response.plot2(
    models  = ProLau_rf,
    Data = get_formal_data(ProLau_models,'expl.var'), 
    show.variables= get_formal_data(ProLau_models,'expl.var.names'),
    do.bivariate = FALSE,
    fixed.var.metric = 'median',
    legend = FALSE,
    display_title = FALSE,
    data_species = get_formal_data(ProLau_models,'resp.var')
  )

gam_eval_strip <- 
  biomod2::response.plot2(
  models  = ProLau_gam,
  Data = get_formal_data(ProLau_models,'expl.var'), 
  show.variables= get_formal_data(ProLau_models,'expl.var.names'),
  do.bivariate = FALSE,
  fixed.var.metric = 'median',
  legend = FALSE,
  display_title = FALSE,
  data_species = get_formal_data(ProLau_models,'resp.var')
)

8.3.1.1.5 ensemble models

## run the ensemble models ----
ProLau_ensemble_models <- 
  BIOMOD_EnsembleModeling(
    modeling.output = ProLau_models,
    em.by = 'all',
    eval.metric = 'TSS',
    eval.metric.quality.threshold = 0.8,
    models.eval.meth = c('TSS','ROC'),
    prob.mean = FALSE,
    prob.cv = TRUE, 
    committee.averaging = TRUE,
    prob.mean.weight = TRUE,
    VarImport = 0
  )

## asses ensemble models quality ----
(ProLau_ensemble_models_scores <- get_evaluations(ProLau_ensemble_models))

8.3.1.1.6 模型投影

## current projections
ProLau_models_proj_current <-  ##各模型分别投影
  BIOMOD_Projection(
    modeling.output = ProLau_models,
    new.env = bioclim_ZA_sub,
    proj.name = "current",
    binary.meth = "TSS",
    output.format = ".img",
    do.stack = FALSE
  )

ProLau_ensemble_models_proj_current <- ##整合模型投影
  BIOMOD_EnsembleForecasting(
    EM.output = ProLau_ensemble_models,
    projection.output = ProLau_models_proj_current,
    binary.meth = "TSS",
    output.format = ".img",
    do.stack = FALSE
  )

## future projections
## load 2050 bioclim variables
bioclim_ZA_2050_BC45 <-
  stack(
    c(
      bio_5  = '../data/worldclim_ZA/worldclim_ZA_2050_BC45_bio_5.asc',
      bio_7  = '../data/worldclim_ZA/worldclim_ZA_2050_BC45_bio_7.asc',
      bio_11 = '../data/worldclim_ZA/worldclim_ZA_2050_BC45_bio_11.asc',
      bio_19 = '../data/worldclim_ZA/worldclim_ZA_2050_BC45_bio_19.asc'
    )
  )
ProLau_models_proj_2050_BC45 <- 
  BIOMOD_Projection(
    modeling.output = ProLau_models,
    new.env = bioclim_ZA_2050_BC45,
    proj.name = "2050_BC45",
    binary.meth = "TSS",
    output.format = ".img",
    do.stack = FALSE
  )
ProLau_ensemble_models_proj_2050_BC45 <- 
  BIOMOD_EnsembleForecasting(
    EM.output = ProLau_ensemble_models,
    projection.output = ProLau_models_proj_2050_BC45,
    binary.meth = "TSS",
    output.format = ".img",
    do.stack = FALSE
  )
## load 2070 bioclim variables
bioclim_ZA_2070_BC45 <-
  stack(
    c(
      bio_5  = '../data/worldclim_ZA/worldclim_ZA_2070_BC45_bio_5.asc',
      bio_7  = '../data/worldclim_ZA/worldclim_ZA_2070_BC45_bio_7.asc',
      bio_11 = '../data/worldclim_ZA/worldclim_ZA_2070_BC45_bio_11.asc',
      bio_19 = '../data/worldclim_ZA/worldclim_ZA_2070_BC45_bio_19.asc'
    )
  )
ProLau_models_proj_2070_BC45 <- 
  BIOMOD_Projection(
    modeling.output = ProLau_models,
    new.env = bioclim_ZA_2070_BC45,
    proj.name = "2070_BC45",
    binary.meth = "TSS",
    output.format = ".img",
    do.stack = FALSE
  )
ProLau_ensemble_models_proj_2070_BC45 <- 
  BIOMOD_EnsembleForecasting(
    EM.output = ProLau_ensemble_models,
    projection.output = ProLau_models_proj_2070_BC45,
    binary.meth = "TSS",
    output.format = ".img",
    do.stack = FALSE
  )

## check how projections looks like ##投影可视化
plot(ProLau_ensemble_models_proj_2070_BC45, str.grep = "EMca|EMwmean")

8.3.1.1.7 计算投影范围变化（加、不变、减）

## compute Species Range Change (SRC) ----
## load binary projections
ProLau_bin_proj_current <- 
  stack( 
    c(
      ca = "Protea.laurifolia/proj_current/individual_projections/Protea.laurifolia_EMcaByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img",
      wm = "Protea.laurifolia/proj_current/individual_projections/Protea.laurifolia_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
    )
  )

ProLau_bin_proj_2050_BC45 <- 
  stack( 
    c(
      ca = "Protea.laurifolia/proj_2050_BC45/individual_projections/Protea.laurifolia_EMcaByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img",
      wm = "Protea.laurifolia/proj_2050_BC45/individual_projections/Protea.laurifolia_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
    )
  )

ProLau_bin_proj_2070_BC45 <- 
  stack( 
    c(
      ca = "Protea.laurifolia/proj_2070_BC45/individual_projections/Protea.laurifolia_EMcaByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img",
      wm = "Protea.laurifolia/proj_2070_BC45/individual_projections/Protea.laurifolia_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
    )
  )

## SRC current -> 2050
SRC_current_2050_BC45 <- 
  BIOMOD_RangeSize(
    ProLau_bin_proj_current,
    ProLau_bin_proj_2050_BC45
  )

SRC_current_2050_BC45$Compt.By.Models

## SRC current -> 2070
SRC_current_2070_BC45 <- 
  BIOMOD_RangeSize(
    ProLau_bin_proj_current,
    ProLau_bin_proj_2070_BC45
  )

SRC_current_2070_BC45$Compt.By.Models

ProLau_src_map <- 
  stack(
    SRC_current_2050_BC45$Diff.By.Pixel, 
    SRC_current_2070_BC45$Diff.By.Pixel
  )
names(ProLau_src_map) <- c("ca cur-2050", "wm cur-2050", "ca cur-2070", "wm cur-2070")

my.at <- seq(-2.5, 1.5, 1)
myColorkey <- 
  list(
    at = my.at, ## where the colors change
    labels = 
      list(
        labels = c("lost", "pres", "abs","gain"), ## labels
        at = my.at[-1] - 0.5 ## where to print labels
      )
  )

rasterVis::levelplot(   ##可视化范围变化；
  ProLau_src_map, 
  main = "Protea laurifolia range change",
  colorkey = myColorkey,
  col.regions=c('#f03b20', '#99d8c9', '#f0f0f0', '#2ca25f'),
  layout = c(2,2)
)

## compute the stratified density of probabilities on SRC ----
##计算不同条件变化下面积变化的概率密度函数；
## the reference projetion
ref <- subset(ProLau_bin_proj_current, "ca")
## define the facets we want to study
mods <- c("GLM", "GBM", "RF", "GAM")
data_set <- c("PA1", "PA2")
cv_run <- c("RUN1", "RUN2", "Full")
## construct combination of all facets
groups <- 
  as.matrix(
    expand.grid(
      models = mods, 
      data_set = data_set, 
      cv_run = cv_run,
      stringsAsFactors = FALSE)
  )
## load all projections we have produced
all_bin_proj_files <- 
  list.files( 
    path = "Protea.laurifolia",  
    pattern = "_TSSbin.img$",
    full.names = TRUE, 
    recursive = TRUE
  )
## current versus 2070 (removed the projections for current and 2050)
current_and_2070_proj_files <- grep(all_bin_proj_files, pattern="2070", value=T)
## keep only projections that match with our selected facets groups
selected_bin_proj_files <- 
  apply(
  groups, 1, 
    function(x){
      proj_file <- NA
      match_tab <- sapply(x, grepl, current_and_2070_proj_files)
      match_id <- which(apply(match_tab, 1, all))
      if(length(match_id)) proj_file <- current_and_2070_proj_files[match_id]
      proj_file
    }
  )
## remove no-matching groups
to_remove <- which(is.na(selected_bin_proj_files))
if(length(to_remove)){
  groups <- groups[-to_remove,]
  selected_bin_proj_files <- selected_bin_proj_files[-to_remove]
}
## build stack of selected projections
proj_groups <- stack(selected_bin_proj_files)
ProbDensFunc(
  initial = as.vector(ref),
  projections = raster::as.matrix(proj_groups),
  groups = t(groups),
  plothist = TRUE,
  resolution = 10,
  cvsn = FALSE
)

8.3.1.1.8 评估建模中不同参数变化对建模结果的影响

library(biomod2)
library(raster)
library(reshape2)
library(ggplot2)
setwd("workdir")
col_vec = c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c'
            ,'#fdbf6f','#ff7f00','#cab2d6','#6a3d9a','#ffff99','#b15928')
col_fun = colorRampPalette(col_vec)
###################################################################################
## PART 1.a
###################################################################################
## get models ##这里加载的是单物种建模的结果；
ProLau_models <- get(load("Protea.laurifolia/Protea.laurifolia.demo1.models.out"))
## get models evaluation scores
ProLau_models_scores <- get_evaluations(ProLau_models)
## rearrange data
ModEval = melt(ProLau_models_scores)
ModEval = ModEval[which(ModEval$Var2 == "Testing.data"), ]
head(ModEval)
##绘制的是不同abseces数据集和不同模型的tss、auc和kappa值差异
ggplot(ModEval, aes(x = Var1, y = value, color = interaction(Var5,Var3))) +
  scale_color_manual("", values = rep(col_fun(4), each = 3)) +
  geom_boxplot(size = 0.8) +
  labs(x = "", y = "") +
  theme(
    axis.text.x = element_text(angle = 90),
    panel.border = element_rect(size = 1, fill = NA),
    panel.grid.major = element_line(linetype = 2, color = "grey", size = 0.8)
  )

###################################################################################
## PART 1.b
###################################################################################
## get ensemble model
ProLau_ensemble_models <- get(load("Protea.laurifolia/Protea.laurifolia.demo1ensemble.models.out"))
## get ensemble model evaluation scores
ProLau_ensemble_models_scores <- get_evaluations(ProLau_ensemble_models)
## rearrange data
EnsModEval = melt(ProLau_ensemble_models_scores)
EnsModEval = EnsModEval[which(EnsModEval$Var2 == "Testing.data"), ]
head(EnsModEval)
##绘制ensemble的整合过程（求mean）中tss和auc变化柱状图
ggplot(EnsModEval, aes(x = Var1, y = value, color = Var1)) +
  scale_color_manual("", values = rep(col_fun(4), each = 3)) +
  geom_boxplot(size = 0.8) +
  labs(x = "", y = "") +
  theme(
    axis.text.x = element_text(angle = 90),
    panel.border = element_rect(size = 1, fill = NA),
    panel.grid.major = element_line(linetype = 2, color = "grey", size = 0.8)
  )

###################################################################################
## PART 2.a
###################################################################################
## get models
ProLau_models <- get(load("Protea.laurifolia/Protea.laurifolia.demo1.models.out"))
## get models variable importance
ProLau_models_var_import <- get_variables_importance(ProLau_models)
## rearrange data
VarImport = melt(ProLau_models_var_import)
head(VarImport)
##绘制的是不同abseces数据集和不同模型的四种环境因子值之间的差异
ggplot(VarImport, aes(x = Var1, y = value, color = interaction(Var4,Var2))) +
  scale_color_manual("", values = rep(col_fun(4), each = 3)) +
  geom_boxplot(size = 0.8) +
  labs(x = "", y = "") +
  theme(
    axis.text.x = element_text(angle = 90),
    panel.border = element_rect(size = 1, fill = NA),
    panel.grid.major = element_line(linetype = 2, color = "grey", size = 0.8)
  )

###################################################################################
## PART 2.b
###################################################################################
## get ensemble model
ProLau_ensemble_models <- get(load("Protea.laurifolia/Protea.laurifolia.demo1ensemble.models.out"))
## get ensemble model variable importance
ProLau_ensemble_models_var_import <- get_variables_importance(ProLau_ensemble_models)
## rearrange data
EnsVarImport = melt(ProLau_ensemble_models_var_import)
EnsVarImport$Var3 = sub(".*EM", "", EnsVarImport$Var3)
EnsVarImport$Var3 = sub("_merged.*", "", EnsVarImport$Var3)
head(EnsVarImport)
ggplot(EnsVarImport, aes(x = Var1, y = value, color = Var3)) +
  scale_color_manual("", values = col_fun(3)) +
  geom_boxplot(size = 0.8) +
  labs(x = "", y = "") +
  theme(
    axis.text.x = element_text(angle = 90),
    panel.border = element_rect(size = 1, fill = NA),
    panel.grid.major = element_line(linetype = 2, color = "grey", size = 0.8)
  )

###################################################################################
## PART 3
###################################################################################
## get some outputs maps
files.names = list.files(path = "Protea.laurifolia/proj_current/individual_projections/"
                         , pattern = "_PA1_"
                         , full.names = TRUE)
files.names = files.names[grep(".img$", files.names)]
files.names = files.names[-grep("TSSbin.img$", files.names)]
stk.outputs = stack(files.names)
names(stk.outputs) = sub("proj_current_Protea.laurifolia_", "", names(stk.outputs))

stk.pts = rasterToPoints(stk.outputs)
stk.pts = as.data.frame(stk.pts)
stk.pts = melt(stk.pts, id.vars = c("x", "y"))
##绘制的是不同abseces数据集和不同模型的概率分布之间的差异
ggplot(stk.pts, aes(x = x, y = y, fill = value/1000)) +
  geom_tile() +
  scale_fill_viridis_c("Probability") +
  facet_wrap(~ variable) +
  labs(x = "", y = "") +
  theme(
    panel.border = element_rect(size = 1, fill = NA),
    panel.grid.major = element_line(linetype = 2, color = "grey", size = 0.8)
  )

## get coefficient of variation
ras.cv = raster("Protea.laurifolia/proj_current/individual_projections/Protea.laurifolia_EMcvByTSS_mergedAlgo_mergedRun_mergedData.img")

ras.cv.pts = rasterToPoints(ras.cv)
ras.cv.pts = as.data.frame(ras.cv.pts)
colnames(ras.cv.pts) = c("x", "y", "CV")
## 绘制不同条件组合下的变异系数单个出图
ggplot(ras.cv.pts, aes(x = x, y = y, fill = CV)) +
  geom_tile() +
  scale_fill_viridis_c() +
  labs(x = "", y = "") +
  theme(
    panel.border = element_rect(size = 1, fill = NA),
    panel.grid.major = element_line(linetype = 2, color = "grey", size = 0.8)
  )

8.1.1.2 多物种建模求多样性

多物种建模最后会得到一个alpha多样性评估。多物种建模一般是针对同一属下物种建模，使用相似的环境变量以帮助模型快速分析；因为建模所需要构建的物种较多，biomod2提供并行运算的参数方案用于加速建模；

多物种建模的原理与单物种建模是类似的，除了最后会将多个物种建模的数据进行整合成一张图；此外多物种建模也需要按照正常的流程进行加载物种分布数据、环境数据和设置不同模型的参数以及运行ensemble模型等等用于数据整合分析

## setup environment ----
setwd('D:/HIDATA/testohase1/all_scripts/biomod2/biomod2_video_multi_species_modelling/workdir')
## load the required packages
library(biomod2)
library(raster)
library(rasterVis)
library(gridExtra)
library(reshape2)

## read data ----
## species occurences data
data <- read.csv('../data/Larus_occ.csv', stringsAsFactors = FALSE)
head(data)
table(data$species)
spp_to_model <- unique(data$species)

## curent climatic variables
stk_current <- 
  raster::stack(
    c(
      bio_1 =  "../data/worldclim_EU/worldclim_EU_bio_1.tif",
      bio_12 = "../data/worldclim_EU/worldclim_EU_bio_12.tif",
      bio_8 =  "../data/worldclim_EU/worldclim_EU_bio_8.tif"
    ),
    RAT = FALSE
  )

## 2050 climatic variables
stk_2050_BC_45 <- 
  raster::stack(
    c(
      bio_1 =  "../data/worldclim_EU/worldclim_EU_2050_BC45_bio_1.tif",
      bio_12 = "../data/worldclim_EU/worldclim_EU_2050_BC45_bio_12.tif",
      bio_8 =  "../data/worldclim_EU/worldclim_EU_2050_BC45_bio_8.tif"
    ),
    RAT = FALSE
  )

## 2070 climatic variables
stk_2070_BC_45 <- 
  raster::stack(
    c(
      bio_1 =  "../data/worldclim_EU/worldclim_EU_2070_BC45_bio_1.tif",
      bio_12 = "../data/worldclim_EU/worldclim_EU_2070_BC45_bio_12.tif",
      bio_8 =  "../data/worldclim_EU/worldclim_EU_2070_BC45_bio_8.tif"
    ),
    RAT = FALSE
  )


## build species modelling wrapper ----
biomod2_wrapper <- function(sp){
  cat("/n> species : ", sp)

  ## get occurrences points
  sp_dat <- data[data$species == sp, ]

  ## formating the data
  sp_format <- 
    BIOMOD_FormatingData(
      resp.var = rep(1, nrow(sp_dat)), 
      expl.var = stk_current,
      resp.xy = sp_dat[, c("long", "lat")],
      resp.name = sp,
      PA.strategy = "random", 
      PA.nb.rep = 2, 
      PA.nb.absences = 1000
    )

  ## print formatting summary
  sp_format

  ## save image of input data summary
  if(!exists(sp)) dir.create(sp)
  pdf(paste(sp, "/", sp ,"_data_formated.pdf", sep="" ))
  try(plot(sp_format))
  dev.off()

  ## define models options
  sp_opt <- BIOMOD_ModelingOptions()

  ## model species
  sp_model <- BIOMOD_Modeling( 
    sp_format, 
    models = c('GLM', 'FDA', 'RF'), 
    models.options = sp_opt, 
    NbRunEval = 2, 
    DataSplit = 70, 
    Yweights = NULL, 
    VarImport = 3, 
    models.eval.meth = c('TSS', 'ROC'),
    SaveObj = TRUE,
    rescal.all.models = FALSE,
    do.full.models = FALSE,
    modeling.id = "demo2"
  )

  ## save some graphical outputs
  #### models scores
  pdf(paste0(sp, "/", sp , "_models_scores.pdf"))
  try(gg1 <- models_scores_graph(sp_model, metrics = c("TSS", "ROC"), by = 'models', plot = FALSE))
  try(gg2 <- models_scores_graph(sp_model, metrics = c("TSS", "ROC"), by = 'data_set', plot = FALSE))
  try(gg3 <- models_scores_graph(sp_model, metrics = c("TSS", "ROC"), by = 'cv_run', plot = FALSE))
  try(grid.arrange(gg1, gg2, gg3))
  dev.off()

  ## build ensemble models
  sp_ens_model <- 
    BIOMOD_EnsembleModeling(
      modeling.output = sp_model,
      em.by = 'all',
      eval.metric = 'TSS',
      eval.metric.quality.threshold = 0.4,
      models.eval.meth = c('TSS','ROC'),
      prob.mean = FALSE,
      prob.mean.weight = TRUE,
      VarImport = 0
    )

  ## do projections
  proj_scen <- c("current", "2050_BC_45", "2070_BC_45")

  for(scen in proj_scen){
    cat("/n> projections of ", scen)

    ## single model projections
    sp_proj <- 
      BIOMOD_Projection(
        modeling.output = sp_model,
        new.env = get(paste0("stk_", scen)),
        proj.name = scen,
        selected.models = 'all',
        binary.meth = "TSS",
        filtered.meth = NULL,
        compress = TRUE,
        build.clamping.mask = FALSE,
        do.stack = FALSE,
        output.format = ".img"
      )

    ## ensemble model projections
    sp_ens_proj <- 
      BIOMOD_EnsembleForecasting(
        EM.output = sp_ens_model,
        projection.output = sp_proj,
        binary.meth = "TSS",
        compress = TRUE,
        do.stack = FALSE,
        output.format = ".img"
      )
  }

  return(paste0(sp," modelling completed !"))
}


## launch the spiecies modelling wrapper over species list ----
if(require(snowfall)){ ## parallel computation
  ## start the cluster
  sfInit(parallel = TRUE, cpus = 5) ## here we only require 4 cpus
  sfExportAll()
  sfLibrary(biomod2)
  ## launch our wrapper in parallel
  sf_out <- sfLapply(spp_to_model, biomod2_wrapper)
  ## stop the cluster
  sfStop()
} else { ## sequencial computation
  for (sp in spp_to_model){
    biomod2_wrapper(sp)
  }
  ## or with a lapply function in sequential model
  ## all_species_bm <- lapply(spp_to_model, biomod2_wrapper)
}

## produce alpha-diversity maps ----

## current conditons
### load binary projections
f_em_wmean_bin_current <- 
  paste0(
    spp_to_model,
    "/proj_current/individual_projections/", 
    spp_to_model, "_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
  )

### sum all projections
if(length(f_em_wmean_bin_current) >= 2){
  ## initialisation
  taxo_alpha_div_current <- raster(f_em_wmean_bin_current[1]) 
  for(f in f_em_wmean_bin_current[-1]){
    taxo_alpha_div_current <- taxo_alpha_div_current + raster(f)
  }
}

## 2050 conditons
### load binaries projections
f_em_wmean_bin_2050 <- 
  paste0(
    spp_to_model,
    "/proj_2050_BC_45/individual_projections/", 
    spp_to_model, "_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
  )

### sum all projections
if(length(f_em_wmean_bin_2050) >= 2){
  ## initialisation
  taxo_alpha_div_2050 <- raster(f_em_wmean_bin_2050[1]) 
  for(f in f_em_wmean_bin_2050[-1]){
    taxo_alpha_div_2050 <- taxo_alpha_div_2050 + raster(f)
  }
}

## 2070 conditons
### load binaries projections
f_em_wmean_bin_2070 <- 
  paste0(
    spp_to_model,
    "/proj_2070_BC_45//individual_projections/", 
    spp_to_model, "_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
  )

### sum all projections
if(length(f_em_wmean_bin_2070) >= 2){
  ## initialisation
  taxo_alpha_div_2070 <- raster(f_em_wmean_bin_2070[1]) 
  for(f in f_em_wmean_bin_2070){
    taxo_alpha_div_2070 <- taxo_alpha_div_2070 + raster(f)
  }
}

## plot the alpha-div maps
levelplot(
  stack(
    c(
      current = taxo_alpha_div_current, 
      in_2050 = taxo_alpha_div_2050, 
      in_2070 = taxo_alpha_div_2070 
    )
  ),
  main = expression(paste("Larus ", alpha, "-diversity")),
  par.settings = BuRdTheme
)

8.1.1.3 并行运算

library(parallel)
cl=makeCluster(20)
clusterEvalQ(cl,{library(R2OpenBUGS);library(DMwR)})
clusterExport(cl,c("unimodal_modelling","bimodal_modelling","unimodal_file","bimodal_file"))
parLapply(cl,1:100,function(i){
  setwd(paste0('sim',as.character(i)))
  d=read.csv("simulations.txt")[,1]
  nam=paste0("sim",i)
  m1=unimodal_modelling(d,nam,wd=getwd(),unimodal_file)
  m2=bimodal_modelling(d,nam,wd=getwd(),bimodal_file)
  setwd('..')
})
stopCluster(cl)

8.3.2 Kuenm(script)

8.3.2.1 Kuenm包的安装说明

### R包说明：


# (ONLY IF NEEDED) installing the package 
## devtools
install.packages("devtools")
## kuenm（在4.02版本下不能安装，在3.6版本下可以安装）
devtools::install_github("marlonecobos/kuenm")

8.3.2.2 kuenm特色功能

## 外推风险分析 mop
    如果在创建最终模型时执行了转移，则可以使用kuenm_mmop函数评估外推的风险。该功能通过比较校准区域与一个或多个生态位模型已转移到的一个或多个区域或方案之间的环境值，计算面向流动性的奇偶性（MOP）层（Owens等人，2013）。使用此功能生成的图层有助于可视化，前提是存在严格的外推风险，并且投影区域或场景与校准区域之间的相似度不同。分析的结果将写入特定数据中的每组变量（图1，深绿色区域）。

## 建模后分析
### 跨多个参数设置的模型统计信息：
计算跨模型的模型统计信息（允许来自不同参数设置的结果
### 适应性和适合区域之间的变化：
确定合适区域的变化和预测的适用性（考虑多种气候模型，排放情景和时间）
### 来自不同来源的栅格方差层：
创建地图，显示每个变异源贡献的变异；
### 来自不同来源的方差的分层划分：
对模型的方差进行分层划分，这些方差来自不同的变化来源。考虑的四个潜在变化源是：模型复制，参数设置，GCM和RCP。

8. 3.3 wallace(shiny)

install.packages("wallace")
library(wallace)
run_wallace()   ##运行run后可直接进入到浏览器界面

## 关于此包的功能说明与解释：
## 此包的可视化功能较为好用，数据分析工作能够很好的利用数据平台展开；其中有用的部分是将ENMeval包和可视化物种分布点以及相关引用问下纳入到shiny展示中，鉴于用户可视化人机交互

8.3.4 NicheToolBox(shiny)

## 安装与运行：
devtools::install_github('luismurao/nichetoolbox')
library(nichetoolbox)
run_nichetoolbox()

## nichetoolbox的安装说明：
https://luismurao.github.io/GSoC/gsoc_final_eval.html

8.3.4.1 2029-49功能

## GUI使用：
### AppSettings部分
加载利基层
工作流程
## 数据部分
搜寻GBIF记录
GBIF数据清理
用户数据
## 使用动态地图进行地理探索
显示经度和纬度数据
## 使用动态地图进行数据策划
定义一个“ M”图
### 保存工作流程
工作流程报告
## 生态位空间（进行聚类分析的步骤）
1.从栅格图层中提取利基值
2.生态位探索
2.1生态位趋势
3.生态位聚类
(## 提供三维bios空间下的物种生态空间投影可视化)
4.生态位相关性
6.生态位建模
    6.1最小体积椭球模型
（##类似于mve，可构建背景空间，so beautiful）
    6.2 Bioclim模型
    6.3 MaxEnt模型
    地理空间中的ENM投影
7.物种分布模型的绩效
    7.1部分ROC
    7.2二进制图
## 参考资料

8.3.5 ODMAP(shiny)

## ODMAP包的使用说明：
ODMAP协议有两个主要目的。首先，它为作者提供了一个清单，详细列出了模型构建和分析的关键步骤。其次，它引入了一种标准的文档编制方法，可确保透明度和可重复性，便于同行评审和专家评估模型质量以及进行荟萃分析。

## 引用论文：
https://onlinelibrary.wiley.com/doi/full/10.1111/ecog.04960

## github主页：
https://github.com/ChrKoenig/ODMAP

## shiny可视化平台：
https://odmap.wsl.ch/

构建组件以及标准

写作构建方式(标准)

8.3.6 enmSDM(script)

## enmSDM简介：
这个软件包是对Robert Hijmans的R流行的`dismo`软件包的补充。它包含一套效率函数，用于准备数据，训练和评估物种分布模型和生态位模型以及比较生态位。

## R包说明文档：
https://github.com/adamlilith/enmSdm

##软件功能类型如下：
## 1、分布点自相关清理
elimCellDups：消除栅格每个像元中的重复点
geoFold：生成地理上不同的k倍
geoThin和geoThinApprox：地理上的薄点

## 2、模型训练
trainByCrossValid：用于trainXYZ在交叉验证折叠中实现某些功能的包装器（另请参见summaryByCrossValid）。
trainBrt：增强的回归树（BRT）
trainCrf：条件回归树（CRF）
trainGam：通用加性模型（GAM）
trainGlm：广义线性模型（GLM）
trainGlmDredge：广义线性模型（GLM）
trainLars：最小角度回归模型（LARS）
trainMaxEnt和trainMaxNet：Maxent模型
trainNs：自然花键（NSs）
trainRf：随机森林（RF）

## 3、模型评估：
aucWeighted：AUC（有/无站点权重）
aucMultiWeighted：AUC的多元版本（有/无网站权重）
contBoyce：连续博伊斯指数（有/无网站权重）
contBoyce2x：连续博伊斯指数的“ 2X覆盖率”版本（有/无网站权重）
fpb：Fpb（有/无站点权重）
thresholdWeighted：将连续预测转换为二进制预测的阈值（有/无站点权重）
thresholdStats：基于阈值预测（具有/不具有站点权重）对性能统计数据进行建模
tssWeighted：真实技能统计（TSS）（有/无站点权重）
modelSize：模型对象中的响应值数量

## 4、生态位重叠
compareNiches：利基重叠指标
compareResponse：比较利基模型响应和单个变量
mop：根据Saupe等人的方法，计算面向流动性的奇偶性，即对多元距离的度量。2012。
nicheOverlap：按照Broennimann等人的方法计算生态位重叠。全球生态与生物地理21：481-497
randPointsRespectingSelf：随机分配地理点，同时大约尊重观察到的点之间的空间自相关结构
randPointsRespectingSelfOther2：随机分配两组地理点，同时大约尊重组之间和组内观察到的空间自相关结构
randPointsBatch：通话randPointsRespectingSelf或randPointsRespectingSelfOther2多次
randPointsBatchExtract：从一组栅格中提取环境以获取使用生成的随机点集 randPointsBatch
randPointsBatchSampled：整理使用生成的所有随机点集 randPointsBatch
randPointsBatchNicheOverlap：计算使用生成的随机点集之间的利基重叠 randPointsBatch

## 5、空间自相关
localSpatialCorrForValues：计算与点或栅格关联的值的局部空间自相关的特定于站点的特征距离
spatialCorrForPoints：计算基于成对距离的地理点之间全局空间自相关的度量
spatialCorrForPointsSummary：空间点的特征簇大小（全局自相关距离）
spatialCorrForPointsPlot：绘制基于成对距离的全局空间自相关度量的观测值和零分布
spatialCorrForPointsWeight：基于全局空间自相关的基于成对距离的度量为点分配权重

## 6、基于最小凸多边形的范围区域
mcpFromPolygons：一组多边形和点中的最小凸多边形
areaFromPointsOrPoly：空间多边形或一组点的面积

## 7、十进制和度分秒形式转换：
decimalToDms：将十进制坐标转换为度-分-秒
dmsToDecimal：将度-分-秒坐标转换为小数

8.3.7 zoon(script)

##zoon包为模块化工作，将工作流定义为模块，然后将模块安装到流程集中，在流程集中批量运行：
图书馆（zoon）


library(zoon)
# 运行工作流程，并指定每种类型的一个模块。
work1 <- workflow(occurrence = UKAnophelesPlumbeus,
                  covariate  = UKAir,
                  process    = OneHundredBackground,
                  model      = LogisticRegression,
                  output     = PrintMap)
# Get a list of modules
GetModuleList()
# Get help on a module
ModuleHelp(LogisticRegression)

## R包使用教程：
https://zoonproject.wordpress.com/tutorials/
## R包使用实例：
https://github.com/lifewatch/sdmpredictors

##  参考文献：
* Species Distribution Modelling has many related names with subtle differences and many similarities such as Ecological or Environmental Niche Modelling, Bio-climate Modelling, Habitat Suitability Modelling, Resource Selection Modelling, Climate Envelope Modelling, Climate Affinity Modelling…

** MaxEnt, R, BioMod, OpenModeller, ModEco, GARP, BioMapper, Canoco,Winbugs, OpenBugs, Domain, ANN, AquaMaps, BioClim, BRT, CSM, CTA, ENFA, Envelope Score, Env Distance, GA, GAM, GBM, Dismo, Glm, GLS, Mahalanobis Distance, MARS, Random Forests, Species, HyperNiche, SRE,SVM , Graf, INLA, BayesComm…

8.3.8 ENMTML(script)

## 关于ENMTML 的解读与说明；
## 参考网址：https://rdrr.io/github/andrefaa/ENMTML/f/README.md
## 论文参考链接：
ENMTML: An R package for a straightforward construction of complex ecological niche models
## 该包是构建将若干篇文献中的方法构建结合起来，然后汇总到一个函数式中；
## 所有的流程均可调参；
## 可同时对多个物种进行批量建模；可并行运行；
### R代码实现：
ENMTML(pred_dir, 
       proj_dir = NULL, 
       result_dir = NULL,
       occ_file, 
       sp, 
       x, 
       y, 
       min_occ = 10,
       thin_occ = NULL, 
       eval_occ = NULL, 
       colin_var = NULL,
       imp_var = FALSE, 
       sp_accessible_area = NULL, 
       pseudoabs_method,
       pres_abs_ratio = 1, 
       part, save_part = FALSE, 
       save_final = TRUE,
       algorithm, 
       thr, 
       msdm = NULL, 
       ensemble = NULL,
       extrapolation = FALSE, 
       cores = 1)
  ## 详细参数说明：
  pred_dir：字符。具有预测变量的目录路径（支持的文件格式为ASC，BILL，TIFF或TXT）
proj_dir：字符。目录路径，其中包含用于项目模型的具有不同区域或时间段的预测变量的文件夹（支持的文件格式为ASC，BILL，TIFF或TXT）。
result_dir：字符。包含要在其中记录模型结果的文件夹的目录路径：
NULL：结果将记录在默认的Result文件夹中，与pred_dir文件夹处于同一级别（用法 result_dir=NULL）。
简单名称：将在与pred_dir文件夹相同的级别上创建具有指定名称的文件夹（例如用法 result_dir="MyFolderName"）
完整路径：将在指定路径（例如result_dir="C:/Users/mypc/Documents/MyFolderName"）创建一个文件夹 。
occ_file：字符。制表符分隔的TXT文件的目录路径，该文件至少包含三列，其中包含有关物种名称以及物种发生的纬度和经度的信息。
sp：字符。列名称，其中包含有关物种名称的信息。
x：字符。列的名称，其中包含有关经度的信息。
y：字符。包含有关纬度信息的列的名称。
min_occ：整数。最小不重复出现次数（少于该次数的物种将被排除）。
thin_occ：字符。对存在的内容执行空间过滤（Thinning，基于spThin包）。对于此扩充，有必要提供一个向量，其中的元素需要具有名称“方法”或“方法”和“距离”（以下更多信息）。三种稀疏方法可用（默认NULL）：
MORAN距离，由Moran变异函数定义。用法 thin_occ=c(method='MORAN')。
像元大小由2倍像元大小定义的距离（Haversine变换）。用法 thin_occ=c(method='CELLSIZE')。
用户定义的用户定义的距离。对于此选项，必须提供一个带有两个值的向量。用法 thin_occ=c(method='USER-DEFINED', ditance='300')。第二个数值是指将用于细化的距离（以km为单位）。因此，distance = 300表示将删除半径300 km之内的所有记录。
eval_occ：字符。制表符分隔的TXT文件的目录路径，带有种类名称，纬度和经度，这三列的列名称必须与occ_file参数中使用的数据库相同。该外部事件数据库将用于外部模型验证（即，它将不用于模型拟合）。（默认NULL）。
colin_var：字符。减少可变共线性的方法：
PCA：对预测变量执行主成分分析，并将主成分用作环境变量。用法 colin_var=c(method='PCA')。
VIF：方差膨胀因子。用法colin_var=c(method='VIF')。
PEARSON：通过Pearson相关选择变量，最大相关阈值必须由用户指定。用法colin_var=c(method='PEARSON', threshold='0.7')。
imp_var：逻辑。对所选算法执行变量重要性和数据以进行曲线响应？（默认为FALSE）
sp_accessible_area：字符。限制每个物种的可访问区域，即用于模型拟合的区域。必须为此参数提供一个向量。实施了三种方法
用于对拟合进行建模的BUFFER区域，该区域由缓冲区界定，缓冲区的宽度等于每种物种在事件对之间的最大距离。用法sp_accessible_area=c(method='BUFFER', type='1')。
用于缓冲拟合的BUFFER区域，该区域由用户定义的宽度大小（以km为单位）缓冲。请注意，此缓冲区的宽度大小将用于所有种类。用法sp_accessible_area=c(method='BUFFER', type='2', width='300')。
屏蔽此方法在于根据物种发生的多边形来划定用于进行模型拟合的区域。例如，有可能基于ecoregion shapefile划定校准区域。对于此选项，有必要告知将用作掩码的文件的路径。可以加载下一个文件格式“ .bil”，“。asc”，“。tif”，“。shp”和“ .txt”。用法sp_accessible_area=c(method='MASK', filepath='C:/Users/mycomputer/ecoregion/olson.shp')..
pseudoabs_method：字符。伪缺席分配方法。必须为此参数提供一个向量。只能选择一种方法。实现以下方法：
RND：在整个用于模型拟合的区域中随机分配的伪缺席。用法pseudoabs_method=c(method='RND')。
ENV_CONST：伪缺失在环境上被限制在Bioclim模型预测的适用性值较低的区域。用法pseudoabs_method = c（method ='ENV_CONST'）。用法pseudoabs_method=c(method='ENV_CONST')。
GEO_CONST：伪缺席的分配是根据地理缓冲区远离发生的事件。对于此方法，有必要提供一个以km为单位表示缓冲区宽度的第二个值。用法pseudoabs_method=c(method='GEO_CONST', width='50')。
GEO_ENV_CONST：伪缺失在环境上受到限制（基于Bioclim模型），但地理分布与基于地理缓冲区的发生位置相距较远。对于此方法，必须提供第二个值，以km为单位表示缓冲区宽度。用法 pseudoabs_method=c(method='GEO_ENV_CONST', width='50')
GEO_ENV_KM_CONST：伪缺席限制在三级过程中；它与GEO_ENV_CONST相似，但有一个额外的步骤，即使用k均值聚类分析在环境空间中分配伪缺失。对于此方法，必须提供第二个值，以km为单位表示缓冲区宽度。用法pseudoabs_method=c(method='GEO_ENV_KM_CONST', width='50')。
pres_abs_ratio：数字。在场与缺席比率（值介于0和1之间）。
部分：性格。用于模型验证的分区方法。只能选择一种方法。必须为此参数提供一个向量。实现以下方法：
引导：随机引导分区。用法part=c(method='BOOT', replicates='2', proportion='0.7')。replicate指重复数，它假定值> = 1。proportion是指用于模型拟合的出现比例，它假定值> 0和<= 1。在此示例中，proportion='0.7'意味着70％的数据将用于模型训练，而30％的数据将用于模型测试。
KFOLD： k折交叉验证中的随机分区。用法part=c(method= 'KFOLD', folds='5')。folds指数据分区的折叠数，它假定值> = 1。
条带：地理分区，其结构是按纬度方式（类型1）或纵向方式（类型2）排列的条带。用法part=c(method= 'BANDS', type='1')。type指乐队的性格。
块：构造为棋盘格的地理分区（又名块交叉验证）。用法part=c(method= 'BLOCK')。
save_part：逻辑。如果为TRUE，函数将保存部分模型的.tif文件，即每个出现的分区创建的模型。（默认FALSE）。
save_final：逻辑。如果为TRUE，则函数将保存最终模型的.tif文件，即包含所有出现数据的文件。（默认TRUE）
算法：字符。构建生态位模型的算法（可以使用多种方法）：
生物： Bioclim
MAH： Mahalanobis
DOM：域
ENF：生态位因子分析
MXS： Maxent简单（仅线性和二次特征，基于MaxNet软件包）
MXD： Maxent的默认功能（所有功能均基于MaxNet软件包）
SVM：支持向量机
GLM：广义线性模型
GAM：广义加法模型
BRT：增强回归树
RDF：随机森林
MLK：最大可能性
GAU：高斯流程使用算法= c（'BIO'，'SVM'，'GLM'，'GAM'，'GAU'）。
thr：字符。用于存在缺席预测的阈值。可以使用多个阈值类型。必须为此参数提供一个向量：
LPT：无遗漏的最高阈值。用法thr=c(type='LPT')。
MAX_TSS：灵敏度和特异性之和最高的阈值。用法。 thr=c(type='MAX_TSS')
MAX_KAPPA： kappa最高的阈值（“ max kappa”）。用法thr=c(type='MAX_KAPPA')。
灵敏度：用户指定的阈值。用法thr=c(type='SENSITIVITY', sens='0.6')。“感觉”是指将使用此适用性值对模型进行二值化。请注意，此方法对所有算法和种类均假设为“敏感”值。
JACCARD： Jaccard最高的阈值。用法thr=c(type='JACCARD')。
SORENSEN： Sorensen最高的阈值。用法thr=c(type='SORENSEN')。
在使用多个阈值类型的情况下，有必要将阈值类型的名称连接起来，例如thr=c(type=c('LPT', 'MAX_TSS', 'JACCARD'))。当SENSITIVITY阈值与其他阈值结合使用时，有必要指定所需的灵敏度值，例如thr=c(type=c('LPT', 'MAX_TSS', 'SENSITIVITY'), sens='0.8')。

msdm：字符。包括模型投影的空间限制。这些方法限制了生态位模型，以减少潜在的预测并使模型更接近物种分布模型。它们分为“先验”和“后验”方法。第一种方法包括将地理层作为模型拟合的预测因子的方法，而后验约束则基于出现和适合性模式来约束模型。此参数仅用一种方法填充，在使用MCP-B方法的情况下，msdm以下面的另一种方式填充（默认NULL）：
a Priori方法（在模型拟合时将图层创建的区域添加为预测变量）：+ XY：创建两个图层的纬度和经度图层。用法msdm=c(method='XY')。+ MIN：创建一个图层，其中包含从每个像元到最近出现的点的距离信息。用法msdm=c(method='MIN')。+ CML：创建一个图层，其中包含从每个像元到所有事件的总距离信息。用法msdm=c(method='CML')。+ KER：在事件数据上创建一个具有高斯核的图层。用法msdm=c(method='KER')。

后验方法：+ OBR：基于出现的限制，使用点之间的距离来排除更合适的补丁（Mendes等人，准备中）。用法msdm=c(method='OBR')。+ LR：下位数，选择最近的25％色标（Mendes等，准备中）。用法msdm=c(method='LR')。+ PRES：仅选择已确认发生数据的补丁（Mendes等，准备中）。用法msdm=c(method='PRES')。+ MCP：排除根据事件数据构建的最小凸多边形（MCP）之外的合适像元。用法msdm=c(method='MCP')。+ MCP-B：在MCP周围创建一个缓冲区（宽度由用户定义，以公里为单位）。用法msdm=c(method='MCP-B', width=100)。在这种情况下width=100 表示将在MCP周围创建一个宽度为100 km的缓冲区。

合奏：性格。用于集成不同算法的方法。可以使用多种方法。必须为此参数提供一个向量。对于SUP，W_MEAN或PCA_SUP方法，有必要为整体参数（例如AUC，Kappa，TSS，Jaccard，Sorensen或Fpb）提供评估指标。（默认为NULL）：
均值：不同模型的简单平均值。用法ensemble=c(method='MEAN')。
W_MEAN：基于模型性能的加权平均值。必须提供评估指标。用法ensemble=c(method='W_MEAN', metric='TSS')。
SUP：最佳模型的平均值（例如，TSS超过平均值）。必须提供评估指标。用法ensemble=c(method='SUP', metric='TSS')。
PCA：执行主成分分析（PCA）并返回第一个轴。用法ensemble=c(method='PCA')。
PCA_SUP：最佳模型的PCA（例如，TSS超过平均值）。必须提供评估指标。用法ensemble=c(method='PCA_SUP', metric='Fpb')。
PCA_THR： PCA仅对适应性值高于所选阈值的那些单元执行。用法ensemble=c(method='PCA_THR')。
如果使用多个集成方法，则必须在参数中连接集成方法的名称，例如ensemble=c(type=c('MEAN', 'PCA'))，ensemble=c(method=c('MEAN, 'W_MEAN', 'PCA_SUP'), metric='Fpb')。

外推逻辑。如果为TRUE，则该函数将基于针对当前条件的面向移动的奇偶性分析（MOP）计算外推。如果使用参数proj_dir，则还将计算其他区域或时间段的外推层（默认值FALSE）。

核心数字。定义要并行运行建模过程的CPU内核数（默认1）。

## 结果文件：
在您的环境变量文件夹之上的一级将创建一个Result文件夹，在内部您将找到一个用于算法结果，整体结果（如果已选择），外推区域的投影和地图的文件夹。还有一些.txt文件（请注意，将在ceratin建模设置下创建一些txt，例如Thresholds_Ensemble.txt）： Evaluation_Table.txt包含模型评估的结果 InfoModeling.txt建模参数的信息 Number_Unique_Occurrences.txt唯一编号每个物种的 出现Occurrences_Cleaned.txt通过在每个网格单元中选择一个出现来过滤出现之后产生的数据集 Occurrences_Filtered.txt通过选择校正后的采样空间偏差（稀疏事件）来过滤事件发生后产生的数据集 Validation_Partition.txt部分模型评估的信息（例如，将模型投影到剩余的30％进行测试时） Thresholds_Algorithm.txt关于阈值的信息用于为每种算法创建在场图（从完整模型的阈值创建在场图） Thresholds_Ensemble.txt有关用于为集合模型创建在场图的阈值的信息

8.3.9 SDMbench(shiny)

## 关于SDMbench的解读：
## git网址：https://github.com/boyanangelov/sdmbench
# 此包的作用是利用shiny构建环境处理体系，比warraace差很多，但是模型的提供的很丰富，不知道是否支持多物种建模；
## 模型是一个整合包，可以用于选择最优模型并利用最优模型进行投影，重要是是可以提供所有方法的参考文献，这就可以降低写作参考文献的难度；

## R实现：
devtools::install_github("boyanangelov/sdmbench")
library(sdmbench)
run_sdmbench()

## 除了gui界面以外，上面的链接中还提供相应的R代码脚本命令：
详细参见，数据的构建方法与传统建模方法一致，但运算workflow比ESDMselsect更简单；

8.3.10 SDMS (script)/Jncc-SDMS

## SDMS的解释与说明:
该包是通过物种继承建模的方式进行，分为两种流式：与biomod2功能类似：可以实现单物种建模和多物种联合建模；
具体脚本参见：收到sdms/vignettes；

## SDMS包的来源说明：
https://github.com/jncc/sdms/tree/master/vignettes

8.3.11 eSDM(shiny)

## eSDM为shiny开发的有gui界面的操作包，可旨在允许用户根据物种分布模型（SDM）在不同的空间比例或使用不同的预测单位创建预测的集合；
## 软件全称：（#原来还是enseable，无趣！）
Ensemble tool for predictions from Species Distribution Modelsgraphical user interface

## 安装R包
devtools::install_github('smwoodman/eSDM',build_vignettes  =  TRUE)
加载R包：
eSDM::eSDM_GUI()

8.3.11 seegSDM(script)

## seegSDM 的解读：
本包的重要优势是可以利用单一model很好的评估该模型建模的优劣以及继续可视化关系；此外，模型还可以构建多次迭代关系，进而确定不确定性，实际上模型的不确定性分布和地理分布关系在环境差异上可能存在某种“定性”的关系，这种关系的发现对于理解物种分布模型的投影以及空间自相关和自下而上的构造理解可能都会起到良好的促进作用；

## 本包的特色功能：
## 可视化模型评价结果，从多个参数值中，如：kappa、auc、sens、spec和pcc等。当然现在更任课的boyce那一套评价方式；
## 计算模型迭代结果的不确定性；
## 教程案例参见：
https://rawgit.com/SEEG-Oxford/seegSDM/master/seegSDM_tutorial.html

8.3.12 MODEL-R(SHINY)

## modleR 简介说明：
modleR是一种工作流，旨在在执行生态位模型（ENM）时自动执行并记录一些常见步骤。给定事件记录和一组环境预测器，setup_sdmdata()可以通过清理重复项，删除没有环境信息的事件并应用某些地理和环境过滤器来准备数据。它还使用交叉验证或引导过程将数据划分为测试集和训练集。do_any()或do_many()使用多种算法拟合生态位模型，其中一些算法已在dismo软件包中实现（Hijmans等，2017），而其他算法则来自R环境中的其他软件包，例如glm，Support Vector Machines（kernlab和e1071））和随机森林（randomForest）。中提供了一种以几种方式连接各个分区的功能final_model()。最后，ensemble_model()根据不同算法组装模型并提供摘要栅格。

## github主页：
https://github.com/Model-R/modleR

## model-R的关键特色：
1、该R包构建了一整套系统流程，主要工作在于数据清理，参数评估（分割参数和重新），构建模型评估流程，基于评估后模型构建ensemble model（##这里构建ensemb的方法采用平均值和分割法构建两种）

2、此R包构建过程中采用的方式和KNENM是类似的，采用文件夹树的形式来构建；

## 关于该包的关联：
## 和rangeModelMetadata有关联：    
https://github.com/cmerow/rangeModelMetadata

8.3.13 ntbox(shiny)

## ntbox是nichetoolbox(见8.3.4)的稳定版，包含nichetoolsbox的所有功能，此外还包含有几项更新后的功能：

### 模型外推分析和不确定性分析：
1、面向移动的奇偶校验（MOP）（Owens等人，2013）
2、多元环境相似表面（MESS）（Elith，Kearney和Phillips 2010）
3、外推检测工具（ExDet）（Mesgaran，Cousens和Webber 2014）